!********************************************************************************************
!       THIS PROGRAM ACCOMPANIES Fried et al (2024)
!
!       "Understanding the Inequality and Welfare Impacts of Carbon Tax Policies"
!
!   This code is used to solve the steady state with the carbon tax
!   The user defines how the carbon tax revenue is used
!   If you use this program, I would kindly ask that you cite my paper.
!*********************************************************************************************

PROGRAM hetero_ss
USE functions_mod
USE grid_interpolation
USE lambda0_optimizer_mod
USE retired_valuer_mod
USE working_valuer_mod
!USE working_valuer_test_mod
USE parameters_mod
USE parameters_energy_mod
USE NEQNF_INT
USE cev_finder_optimizer_mod

IMPLICIT NONE
INTERFACE
    function randu(iseed)
        INTEGER, intent(inout)   ::iseed
    DOUBLE PRECISION   ::randu
end function randu
END INTERFACE
include 'mpif.h'

parameter send_data_tag = 2001, return_data_tag = 2002

INTEGER            :: con_pos,kstart, capital_positive, tax_loop, non_clear,temp,kit,i,s,it,it1,it4,it_tax
INTEGER            :: ac,out_loop,e_it,rebate_loop
INTEGER            :: sim_shocks(J,simulationsN),sim_types(J,simulationsN)
INTEGER            :: taxes_clear(max_outer_loops+1),krow
INTEGER            :: shocks_start1,shocks_start2,routine
INTEGER            :: num_procs_run,seed,it2,it3,extra,optimal_baseline

INTEGER            :: my_id, root_process, ierr, status(MPI_STATUS_SIZE)
INTEGER            :: num_procs, an_id,vary_capital,update_lambda0_ig
INTEGER            :: num_procs_run2,share_it,share_it2,share_it3
INTEGER            :: start_row,end_row,increment,large_inc,increment2,large_inc2
INTEGER            :: bad_combination_inner_nostop,bad_combination_outer,bad_combination_inner_stop


DOUBLE PRECISION   :: b,pricefactor,buf(J),xguess(3),solver_out(3),K1_share,N1_share,rate_use,ypp_temp
DOUBLE PRECISION   :: curv,diff_parameter,diff_parameter_v(6),survive(J), dying(J),types(typesN)
DOUBLE PRECISION   :: growth,  pop_temp1(J),hours_jnr,sim_earnings(J,simulationsN)
DOUBLE PRECISION   :: sim_labor(J,simulationsN),sim_capital(J+1,simulationsN),sim_consumption(J,simulationsN)
DOUBLE PRECISION   :: sim_energy(J,simulationsN),sim_utility(J,simulationsN),&
    sim_lifetime_total_income(simulationsN),sim_lifetime_after_tax_income(simulationsN),&
    sim_lifetime_total_income_nosurvive(simulationsN),sim_lifetime_after_tax_income_nosurvive(simulationsN),&
    sim_rebate(J,simulationsN)

DOUBLE PRECISION   :: agg_labor(J,typesN),agg_capital(J+1,typesN),agg_consumption(J,typesN),agg_energy(J,typesN)
DOUBLE PRECISION   :: kprime(typesN,shocksN,kgridN,J-1),labor(typesN,shocksN,kgridN,Jnr)
DOUBLE PRECISION   :: dummy,sim_v(simulationsN),sim_total_income(J,simulationsN),sim_after_tax_income(J,simulationsN)
DOUBLE PRECISION   :: v(typesN,shocksN,kgridN,J+1),v_retired(typesN,shocksN,kgridN,J+1),sim_consumption_tilda(J,simulationsN)
DOUBLE PRECISION   :: v_in(typesN,shocksN,kgridN),aggregates_start(31),xguesses_check(3),check_dummy,check_dummy_out(3)
DOUBLE PRECISION   :: ctot80(kgridN,1),e80(kgridN,1),c80(kgridN,1), con(kgridN),con_dummy,sim_labor_tax_new(J,simulationsN),sim_labor_tax_old(J,simulationsN)
DOUBLE PRECISION   :: ss(typesN,(shocksN)/2),ss1(typesN,(shocksN)/2),ss_start(typesN,(shocksN)/2)

DOUBLE PRECISION   :: energy_miss,ind_eng_tax,lump_sum_ig,giniWelSS,lifeWelBar,giniadj,&
    giniIncSS_nosurvive,giniIncSS_aftertax_nosurvive,giniIncSS,giniIncSS_aftertax,giniWelSS_start,giniIncSS_start
DOUBLE PRECISION   :: c_upper(kgridN),avg_earnings1,tax_dummy,labor_tax_dummy,capital_tax_dummy
DOUBLE PRECISION   :: gs,increment_real,increment_real2,hours_jnr_start,hours_jnr_end

DOUBLE PRECISION   :: Ce_for_share_start,C_for_share_start,pe_start,sim_utility_c(J,simulationsN),sim_utility_l(J,simulationsN),&
r_start,lambdak_start

DOUBLE PRECISION   :: rebate_slope_index(slopeN),rebate_level_ig,rebate_level1,rebate_level_index(levelN),&
    rebate_slope_l,rebate_slope_h,rebate_level_l,rebate_level_h,&
    progressivity_l,progressivity_h,progressivity_index(progressivityN),lambda1_benchmark

!double precision aggregates(135,levelN*slopeN*progressivityN+progressivityN+1),parameters(slopeN*levelN*progressivityN+1+progressivityN,3)
!double precision aggregates(135,levelN*slopeN*progressivityN+progressivityN+1),parameters(slopeN*levelN*progressivityN,3)

double precision aggregates(135,levelN*slopeN*progressivityN+1),parameters(slopeN*levelN*progressivityN+1,3)

DOUBLE PRECISION   :: sim_total_income_sorted(J,simulationsN),sim_lifetime_after_tax_income_sorted(1,simulationsN)

DOUBLE PRECISION   :: CE_end,C_end,Ce_start,C_start,sim_consumption_sorted_start(J,simulationsN),&
    sim_labor_sorted_start(J,simulationsN),sim_energy_sorted_start(J,simulationsN)

DOUBLE PRECISION   :: cev_level,&
    cev_total_alt,cev_total_collapsed,cev_total






DOUBLE PRECISION   :: sim_labor_dummy(1,simulationsN),sim_labor_dummy_start(1,simulationsN),&
    sim_consumption_dummy(1,simulationsN),sim_consumption_dummy_start(1,simulationsN)


DOUBLE PRECISION   :: perc_better
DOUBLE PRECISION   :: hbar,partprod, K_ig,E,Ec,Ec1,Ep, Ec_ig,N_ig, tr_ig, K, N, N1,N2,K1,K2,Ae,Atfp,sigmaa,gov_spend
DOUBLE PRECISION   :: K1_prime,N1_prime,Y, w, r,kgridl, kgridu,  tr1,sigmav,sigmaeps,rho
DOUBLE PRECISION   :: a,a1,tax_diff,tax_diff_keep,govt_taxes,govt_taxes1,govt_taxes2,govt_taxes3,govt_taxes4,tol,pop_start,tr
DOUBLE PRECISION   :: shocks_temp((shocksN)/2),stat_dist_temp((shocksN)/2),shock_pi_temp((shocksN)/2,(shocksN)/2)
DOUBLE PRECISION   :: shocks(shocksN),shock_pi(shocksN,shocksN,Jnr-1,typesN),shock_cumpi(shocksN,shocksN,J,typesN)
DOUBLE PRECISION   :: avg_earnings,prices(11)
DOUBLE PRECISION   :: random(J,simulationsN)
DOUBLE PRECISION   :: avg_earnings_ss1(typesN,(shocksN)/2),&
    avg_earnings_ss_ig(typesN,(shocksN)/2),avg_earnings_ss(typesN,(shocksN)/2)
DOUBLE PRECISION   :: ss_count(typesN,(shocksN)/2)

DOUBLE PRECISION   :: sim_consumption_start(J,simulationsN),sim_labor_start(J,simulationsN),sim_energy_start(J,simulationsN),&
    sim_utility_start(J,simulationsN),sim_utility_c_start(J,simulationsN),sim_utility_l_start(J,simulationsN),&
    sim_total_consumption_start(J,simulationsN),sim_labor_earnings_start(J,simulationsN),sim_labor_tax_start(J,simulationsN),sim_capital_earnings_start(J,simulationsN),sim_capital_tax_start(J,simulationsN),&
    sim_total_after_tax_earnings_start(J,simulationsN),sim_lifetime_after_tax_income_start(simulationsN)

DOUBLE PRECISION, allocatable, dimension(:,:) :: zeros, population,sim_labor_alloc,sim_capital_alloc,dummy_zeros
DOUBLE PRECISION, allocatable, dimension(:,:) :: sim_consumption_alloc,sim_earnings_alloc,sim_energy_alloc,&
    sim_total_income_alloc,sim_after_tax_income_alloc
DOUBLE PRECISION, allocatable, dimension(:) :: zeros2,sim_v_alloc,share_guesses,share_guesses2,share_guesses3
DOUBLE PRECISION, allocatable, dimension(:,:,:) :: v_alloc,kprime_alloc,labor_alloc


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MPI Parameters

root_process = 0
call MPI_INIT(ierr)
call MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr)
call MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr)
num_procs_run=num_procs
num_procs_run2=num_procs
!Increments for simulations
if (simulationsN<(num_procs-1)) then
    num_procs_run2=simulationsN+1
end if
increment_real2=simulationsN/(num_procs_run2-1.0D0)
increment2=floor(increment_real2)
large_inc2=anint((increment_real2-increment2)*(num_procs_run2-1))

!Increments for value function calculation

if (kgridN<(num_procs-1)) then
    num_procs_run=kgridN+1
end if
increment_real=kgridN/(num_procs_run-1.0D0)
increment=floor(increment_real)
large_inc=anint((increment_real-increment)*(num_procs_run-1))


CALL ERSET(4,1,0)
CALL ERSET(3,1,0)
CALL ERSET(2,1,0)

household=0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Computational Parameters
a=0.04D0
a1=0.04D0
tol=0.00025D0
!tol=0.94D0
allocate(zeros2(max_outer_loops+1))
zeros2=0
taxes_clear=zeros2
deallocate(zeros2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Utility Parameters
sigma1=2.0D0
sigma2=0.5D0

if(household==0) then
    !Individual
    gamma1=0.99074D0
!    ebar=0.001324D0*2.0D0
    ebar=0.001324D0
!    ebar=0.00D0
    beta=0.99502D0
    chi=73.3
    sigma2=0.5D0
    kgridl=-0.157D0
elseif (household==1) then
    !Household
    gamma1=0.99275D0
    ebar=0.0068D0
    beta=1.0085D0
    chi=90.0D0
    sigma2=0.55D0
    kgridl=0.0D0
else
    print*,"Set Household or Individual"
end if


!!!Set this to 1 if using unconstrained optimal as the baseline; 0 otherwise
optimal_baseline=1


rebate_to_old=1
negative_taxes=0
!Can rebate lead to negative effective taxes
rebate_negative_tax=1
!!!! 1 labor only, 2 capital only, 3 all
labor_capital_rebate_base=3


!No carbon tax but rebate revenue
no_carbon_tax=0


!Make rebate level grid capital grid (1=added as rebate level,2=added as rebate slope)
vary_capital=1



no_increase=1
!Doesn't allow labor tax to be larger than baseline
!-1 no longer constrains labor tax to be lower than previous
extra=0
!have the extra combinations on top, make sure decleration is right above


routine=4
!1 lump_sum_rebate
!2 throw_away
!3 cap_tax_rebate
!4 labor_rebate
!5 Rebate level

!Updated lambda0 when new parameters for IG
update_lambda0_ig=1
if(routine.NE.4) then
    update_lambda0_ig=0
end if

!Create tax that is 0 up to progressivity, then old tax after (routine should be 3)
flat_tax=0
if(flat_tax==1 .and. routine .NE. 3) then
print*,"UHOHOHOHOH flat tax on "
end if



!
rebate_level_l=0.071D0
rebate_level_h=0.079D0
rebate_level_index(1)=rebate_level_l
if(levelN>1) then
    rebate_level_index(levelN)=rebate_level_h
    do it =2,levelN-1
        rebate_level_index(it)=rebate_level_index(it-1)+(rebate_level_h-rebate_level_l)/real(levelN-1.0D0)
    end do
end if
!

progressivity_l=0.211D0-0.031D0
progressivity_h=0.219D0-0.031D0
progressivity_index(1)=progressivity_l
progressivity_index(progressivityN)=progressivity_h
if(progressivityN==1) then
    progressivity_index=0.0D0
else
    do it =2,progressivityN-1
        progressivity_index(it)=progressivity_index(it-1)+(progressivity_h-progressivity_l)/real(progressivityN-1.0D0)
    end do
end if

it2=1

!Check flat grid first
if(vary_capital==1) then
    parameters(it2,1)=0.360d0
    parameters(it2,1)=rebate_level_l
    parameters(it2,2)=0.0d0
    parameters(it2,3)=progressivity_l
elseif(vary_capital==2) then
    parameters(it2,2)=0.360d0
    parameters(it2,1)=0.0D0
else
    parameters(it2,1)=0.0D0
    parameters(it2,2)=0.0d0
end if
if(my_id==1) then
    print*,"out",parameters(1,:),progressivity_l,it2,parameters(it2,3)
end if
!if(flat_tax==1) then
!    parameters(it2,3)=progressivity_l
!else
!    parameters(it2,3)=0.0d0
!end if
it2=it2+1
!
!do it3=1,progressivityN
!if(vary_capital==1) then
!    parameters(it2,1)=0.360d0
!    parameters(it2,2)=0.0D0
!elseif(vary_capital==2) then
!    parameters(it2,2)=0.360d0
!else
!    parameters(it2,1)=0.0d0
!    parameters(it2,2)=0.0D0
!end if
!    parameters(it2,3)=progressivity_index(it3)
!    it2=it2+1
!end do

!it3=1
do it3=1,progressivityN
!it=1
!it1=1
    do it=1,levelN


    do it1=1,slopeN
if(my_id==1) then
print*,"it2",it2,parameters(it2,3),parameters(1,3)
end if
        parameters(it2,1)=rebate_level_index(it)
        parameters(it2,2)=0.0D0
        parameters(it2,3)=progressivity_index(it3)
if(my_id==1) then
print*,"it2a",it2,parameters(it2,3),parameters(1,3)
end if
        it2=it2+1
    end do
    end do
end do

if(my_id==1) then
print*,"parms3",parameters(:,3)
print*,"prog_index",progressivity_index
print*,"parms1",parameters(:,1)
end if





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Productivity Parameters
!Individual
hbar=0.00D0
partprod=1.0D0

if(household==0) then
epsilon1(:,1) = (/1.0D0, 1.05037955556458D0, 1.10168574952967D0, 1.15381030928342D0, 1.20663608757434D0, 1.26003734304819D0, 1.3138800980071D0, 1.36802257372016D0, 1.42231570281972D0, 1.47660371749329D0, 1.53072481133652D0, 1.584511871879D0, 1.63779327994293D0, 1.6903937711557D0, 1.74213535412359D0, 1.79283827899664D0, 1.84232204942572D0, 1.89040647024333D0, 1.93691272260084D0, 1.9816644577762D0, 2.02448890043734D0, 2.06521795181501D0, 2.10368928301267D0, 2.13974740856428D0, 2.1732447303485D0, 2.20404254208235D0, 2.23201198484835D0, 2.25703494445744D0, 2.27900488191109D0, 2.29782758879663D0, 2.31342186012357D0, 2.32572007787805D0, 2.33466869942757D0, 2.34022864584003D0, 2.34237558617614D0, 2.34110011486113D0, 2.33640782032644D0, 2.32831924422039D0, 2.31686973160448D0, 2.30210917366396D0, 2.28410164555347D0, 2.26292494305636D0, 2.23867002274592D0, 2.21144035128553D0, 2.18135117038086D0, 2.14852868468884D0/)
elseif(household==1) then
epsilon1(:,1) = (/1.0D0, 1.11972993442023D0, 1.23163188204661D0, 1.33628274802356D0, 1.43423225535758D0, 1.52600294491726D0, 1.61209017543323D0, 1.69296212349824D0, 1.76905978356708D0, 1.84079696795664D0, 1.90856030684589D0, 1.97270924827585D0, 2.03357605814965D0, 2.09146582023247D0, 2.14665643615159D0, 2.19939862539634D0, 2.24991592531814D0, 2.29840469113051D0, 2.34503409590901D0, 2.3899461305913D0, 2.4332556039771D0, 2.47505014272822D0, 2.51539019136855D0, 2.55430901228404D0, 2.59181268572273D0, 2.62788010979474D0, 2.66246300047226D0, 2.69548589158955D0, 2.72684613484296D0, 2.75641389979091D0, 2.7840321738539D0, 2.8095167623145D0, 2.83265628831737D0, 2.85321219286924D0, 2.87091873483891D0, 2.88548299095726D0, 2.89658485581726D0, 2.90387704187394D0, 2.90698507944442D0, 2.90550731670788D0, 2.89901491970559D0, 2.88705187234091D0, 2.86913497637924D0, 2.84475385144809D0, 2.81337093503703D0, 2.77442148249772D0/)
else
end if
!Firm
alpha1 = 0.3D0
alpha2 = 1.0D0-0.403D0
theta=0.03D0
delta = 0.079062D0
Atfp=1.0D0
Ae=Atfp*1.0D0

!taue=pe*0.326144D0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Heterogeneity Parameters
if(household==0) then
    rho=0.958D0
    sigmaa=0.065D0**0.5D0
    sigmav=(0.017D0/(1.0D0-rho**2.0D0))**0.5D0
    sigmaeps=0.081D0**0.5D0
elseif(household==1) then
    rho=0.972D0
    sigmaa=0.022D0**0.5D0
    sigmav=(0.012D0/(1.0D0-rho**2.0D0))**0.5D0
    sigmaeps=0.093D0**0.5D0
else
    print*,"Set Household"
end if



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Calculate Heterogeneity
do i =1,typesN
    types(i)=exp(-(sigmaa)+(real(i-1)/real(typesN-1))*2.0D0*sigmaa)
end do

! Stochastic parameters
CALL rouwenhorst(rho,sigmav,shock_pi_temp,shocks_temp,stat_dist_temp,(shocksN)/2)


shocks_temp=exp(shocks_temp)
shocks_temp=shocks_temp/(sum(stat_dist_temp*shocks_temp))

shocks(1:int(shocksN/2))=shocks_temp*exp(-sigmaeps)
shocks(int(shocksN/2)+1:int(shocksN))=shocks_temp*exp(sigmaeps)

if(my_id==1) then
print*,"shocks",shocks
end if
shocks_start1=int((((shocksN)/2)+1)/2)
shocks_start2=int(shocks_start1+((shocksN)/2))

!E to E
do it1=1,Jnr-1
do it2=1,typesN
shock_pi(1:int((shocksN)/2),1:int((shocksN)/2),it1,it2)=shock_pi_temp/2.0D0
shock_pi(int((shocksN)/2+1):(shocksN),1:int((shocksN)/2),it1,it2)=shock_pi_temp/2.0D0
shock_pi(1:int((shocksN)/2),int((shocksN)/2+1):(shocksN),it1,it2)=shock_pi_temp/2.0D0
shock_pi(int((shocksN)/2+1):(shocksN),int((shocksN)/2+1):(shocksN),it1,it2)=shock_pi_temp/2.0D0
end do
end do





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Population Parameters
growth = 0.011D0
if(household==0) then
pop_temp1 =(/197316.D0, 197141.D0, 196959.D0, 196770.D0, 196580.D0, 196392.D0, 196205.D0, 196019.D0, 195830.D0, 195634.D0, 195429.D0, 195211.D0, 194982.D0, 194739.D0, 194482.D0, 194211.D0, 193924.D0, 193619.D0, 193294.D0, 192945.D0, 192571.D0, 192169.D0, 191736.D0, 191271.D0, 190774.D0, 190243.D0, 189673.D0, 189060.D0, 188402.D0, 187699.D0, 186944.D0, 186133.D0, 185258.D0, 184313.D0, 183290.D0, 182181.D0, 180976.D0, 179665.D0, 178238.D0, 176689.D0, 175009.D0, 173187.D0, 171214.D0, 169064.D0, 166714.D0, 164147.D0, 161343.D0, 158304.D0, 155048.D0, 151604.D0, 147990.D0, 144189.D0, 140180.D0, 135960.D0, 131532.D0, 126888.D0, 122012.D0, 116888.D0, 111506.D0, 105861.D0, 99957.D0, 93806.D0, 87434.D0, 80882.D0, 74204.D0, 67462.D0, 60721.D0, 54053.D0,	47533.D0,	41241.D0, 35259.D0, 29663.D0, 24522.D0, 19890.D0, 15805.D0, 12284.D0, 9331.D0, 6924.D0, 5016.D0, 3550.D0, 2454.D0 /)
do i=1,(J-1)
   survive(i) = pop_temp1(i+1)/pop_temp1(i)
end do
survive(J) = 0.0D0
adult_equivalence=1.0D0
else
survive=(/1.0D0, 0.999443429520098D0, 0.999453048952094D0, 0.999474759951225D0, 0.999507517658172D0, 0.999550180294631D0, 0.99958688938571D0, 0.999616966794992D0, 0.999637259341625D0, 0.999646435563957D0, 0.999642137784891D0, 0.999638705313936D0, 0.999628264835928D0, 0.999611266702816D0, 0.999587623092204D0, 0.99955347402391D0, 0.999517484919676D0, 0.999474904683543D0, 0.999424216406889D0, 0.999368982231295D0, 0.999311637585569D0, 0.999247262201619D0, 0.99917475362204D0, 0.999096164795857D0, 0.999015366949988D0, 0.99892772761533D0, 0.998832359878012D0, 0.998729478390282D0, 0.998631176268071D0, 0.998535878125485D0, 0.998443264975948D0, 0.998338654631507D0, 0.998221906285934D0, 0.998078357864358D0, 0.997913954710353D0, 0.997725564043014D0, 0.997510605728885D0, 0.997277513805754D0, 0.99702940364917D0, 0.996758720419625D0, 0.996463887451456D0, 0.996136496024502D0, 0.995758248748826D0, 0.995343091935131D0, 0.994874730958277D0, 0.994354016314192D0, 0.993755337681475D0, 0.993065683870379D0, 0.992293047744931D0, 0.991433174915763D0, 0.990457448940283D0, 0.989303075213917D0, 0.987964810901374D0, 0.986475935412401D0, 0.984833077329247D0, 0.982973697702629D0, 0.980769120672249D0, 0.978166137633033D0, 0.975176536129091D0, 0.971754588528722D0, 0.967817497628447D0, 0.963253722024177D0, 0.957955278732916D0, 0.951933505633426D0, 0.945105096525417D0, 0.937431585348977D0, 0.92887640431547D0, 0.919444350437676D0, 0.909124748145856D0, 0.89794551325945D0, 0.885974030114678D0, 0.873304919941224D0, 0.859956433479889D0, 0.842665940409714D0, 0.824214308977887D0, 0.804792055183895D0, 0.785116734289799D0, 0.765756410717599D0, 0.747202819257837D0, 0.729458167443388D0, 0.713669733378468D0/)
survive(J) = 0.0D0
adult_equivalence=(/1.21998398833892D0, 1.26899136580154D0, 1.32056246723579D0, 1.3734138700101D0, 1.42641978111077D0, 1.47860215039009D0, 1.52912105574409D0, 1.57726536021961D0, 1.62244364105098D0, 1.66417539062616D0, 1.70208248938226D0, 1.73588095063074D0, 1.76537293731188D0, 1.79043905067891D0, 1.81103089091152D0, 1.82716388965886D0, 1.83891041451209D0, 1.84639314540638D0, 1.84977872295234D0, 1.84927166869702D0, 1.84510857731432D0, 1.83755258072496D0, 1.8268880841459D0, 1.81341577406915D0, 1.79744789817028D0, 1.77930381714618D0, 1.75930582848249D0, 1.73777526215036D0, 1.71502884823281D0, 1.69137535648057D0, 1.66711250779728D0, 1.64252415765435D0, 1.61787775143512D0, 1.5934220517087D0, 1.5693851374331D0, 1.54597267508807D0, 1.52336646173711D0, 1.50172324001931D0, 1.48117378507042D0, 1.46182226337349D0, 1.44374586353911D0, 1.42699469901482D0, 1.41159198272448D0, 1.39753447363654D0, 1.3847931952624D0, 1.37331442608375D0, 1.36302096190969D0, 1.35381365016324D0, 1.34557319609731D0, 1.33816224094031D0, 1.33142771197095D0, 1.32520344452275D0, 1.31931307591802D0, 1.31357321133145D0, 1.30779686158245D0, 1.30179715285813D0, 1.29539130836514D0, 1.28840490191053D0, 1.28067638341311D0, 1.27206187634362D0, 1.26244024709437D0, 1.25171844627878D0, 1.23983712196005D0, 1.22677650480938D0, 1.21256256519385D0, 1.19727344219358D0, 1.18104614454827D0, 1.16408352353373D0, 1.14666151776743D0, 1.12913666994348D0, 1.11195391549768D0, 1.09896503226485D0, 1.08988645284858D0, 1.0816144615019D0, 1.07430423674565D0, 1.0681188559846D0, 1.06322914430733D0, 1.0598135018135D0, 1.0580577094685D0, 1.05815471348597D0, 1.06030438823754D0/)
!adult_equivalence=1.0D0
end if

dying = 1.D0 - survive
allocate(population(J,2))
population=0.0D0
population(1,1) = 1.0D0
do i=2,(J)
   population(i,1) = survive(i-1)*population((i-1),1)/(1+growth)
end do
pop_start=population(1,1)/sum(population(1:J,1))
population(:,2) = population(:,1)/sum(population(:,1))



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Grids
kgrid(1)=kgridl
kgrid(kgridnegN+1)=0.0D0
curv = 1.5D0
do ac=2,kgridnegN+1
	kgrid(kgridnegN-ac+2)=kgrid(kgridnegN+1)+(kgridl)*((ac-1.0D0)/(kgridnegN))**curv
end do

kgridu = 25.0D0
curv = 1.35D0
do ac=kgridnegN+2,kgridN
	kgrid(ac)=kgrid(kgridnegN+1)+(kgridu-kgrid(kgridnegN+1))*((ac-kgridnegN-1.0D0)/(kgridN-kgridnegN-1.0D0))**curv
end do
kgrid(kgridN)=kgridu
capital_positive = minloc(kgrid,1, mask = kgrid.ge.0)
kstart=capital_positive




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Social Security Payment Information
b=0.45D0

!!!!!!!!!!!!!!!!!!!!!!!!!!!! Government Parameters
gs=.157D0
govt_taxes=0.0D0
lambda1_benchmark=.031D0
lambda1=.031D0

!1 lump_sum_rebate
!2 throw_away
!3 cap_tax_rebate
!4 labor_rebate
if (routine==1 .or. routine == 2 .or. routine==4 .or. routine==5) then
    lambdak=.36D0
!    lambdak=.3125D0

elseif (routine==3 .and. ebar.ne. 0.0D0) then
!    lambdak=0.276263951D0
    lambdak=.36D0
else
    print*,"Not Specified Right"
end if
!lambdak=.0D0

    OPEN(227, FILE="/<filepath starting steady state>/aggregates_output.dat",FORM="formatted")
    read(227,887) (aggregates_start)
    close(227)

    OPEN(227, FILE="/<filepath starting steady state>/ss_output.dat",FORM="formatted")
    read(227,887) (ss_start)
    close(227)

    OPEN(227, FILE="/<filepath starting steady state>/avg_earnings_ss_output.dat",FORM="formatted")
    read(227,887) (avg_earnings_ss_ig)
    close(227)

gov_spend=aggregates_start(7)
lambda0=aggregates_start(6)

Ce_for_share_start=aggregates_start(14)
C_for_share_start=aggregates_start(15)
pe_start=aggregates_start(16)
avg_earnings=aggregates_start(17)
avg_earnings_scale=(avg_earnings)*(1.0D0-aggregates_start(8)/2.0D0)

!!!No tax
!taue=0.0D0

!!Standard tax
taue=pe_start*0.482D0

!Small tax
!taue=pe_start*0.482D0*0.75D0

!!large tax
!taue=pe_start*0.482D0*1.25D0
lambda0_start=aggregates_start(6)
lambda1_start=aggregates_start(13)
!avg_earnings_ss_ig=aggregates_start(28)
avg_earnings_ss=avg_earnings_ss_ig
rebate_level_ig=0.0D0
r_start=aggregates_start(2)
lambdak_start=aggregates_start(12)

if(no_carbon_tax==1) then

print*,"NOT TAXING CARBON"
!Base on optimal
if(optimal_baseline==0) then
gov_spend=aggregates_start(7)-(0.006164949D0+0.025624986D0)*taue
else
gov_spend=aggregates_start(7)-(0.006167261D0+0.026070431D0)*taue
end if
taue=0.0D0
end if

if(my_id==1) then
print*,"here2"
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Initial Guesses
!Capital rebate
!K_ig=1.780835216D0
!N_ig=0.525966934D0
!Ec_ig=0.006223467D0
!tr_ig=0.023550899D0


!!Optimal
K_ig=1.7406D0
N_ig=0.523D0
Ec_ig=0.006167D0
tr_ig=0.0223D0





if(routine==4 .or. routine ==2) then
    lambda0=aggregates_start(6)
end if
!end do
if (optimal_baseline==0) then
 OPEN(227, FILE="/<filepath starting steadystate>/sim_labor_output.dat",FORM="formatted")
    read(227,887) (sim_labor_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steadystate>/sim_energy_output.dat",FORM="formatted")
    read(227,887) (sim_energy_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steadystate>/sim_consumption_output.dat",FORM="formatted")
    read(227,887) (sim_consumption_start)
    close(227)
    close(227)
 OPEN(227, FILE="/<filepath starting steadystate>/sim_labor_earnings_output.dat",FORM="formatted")
    read(227,887) (sim_labor_earnings_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steadystate>/sim_labor_tax_output.dat",FORM="formatted")
    read(227,887) (sim_labor_tax_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steadystate>/sim_capital_earnings_output.dat",FORM="formatted")
    read(227,887) (sim_capital_earnings_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steadystate>/sim_capital_tax_output.dat",FORM="formatted")
    read(227,887) (sim_capital_tax_start)
    close(227)

sim_total_after_tax_earnings_start=sim_labor_earnings_start+sim_capital_earnings_start-sim_capital_tax_start-sim_labor_tax_start
else
 OPEN(227, FILE="/<filepath starting steadystate>/unconstrained_optimal/sim_labor_output.dat",FORM="formatted")
    read(227,887) (sim_labor_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steadystate>/unconstrained_optimal/sim_energy_output.dat",FORM="formatted")
    read(227,887) (sim_energy_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steadystate>/unconstrained_optimal/sim_consumption_output.dat",FORM="formatted")
    read(227,887) (sim_consumption_start)
    close(227)
    close(227)
 OPEN(227, FILE="/<filepath starting steadystate>/unconstrained_optimal/sim_labor_earnings_output.dat",FORM="formatted")
    read(227,887) (sim_labor_earnings_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steadystate>/unconstrained_optimal/sim_labor_tax_output.dat",FORM="formatted")
    read(227,887) (sim_labor_tax_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steadystate>/unconstrained_optimal/sim_capital_earnings_output.dat",FORM="formatted")
    read(227,887) (sim_capital_earnings_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steadystate>/unconstrained_optimal/sim_capital_tax_output.dat",FORM="formatted")
    read(227,887) (sim_capital_tax_start)
    close(227)

sim_total_after_tax_earnings_start=sim_labor_earnings_start+sim_capital_earnings_start-sim_capital_tax_start-sim_labor_tax_start
end if


!change
if(routine==2 .or. routine == 3 .or. routine ==4 .or. routine==5) then
    lump_sum_ig=0.00D0
else if(routine==1) then
    lump_sum_ig=0.05D0
else
    print*,"Not Specified Right"
end if



!Set Up population for simulations
seed=53415691
do s=1,simulationsN
do i=1,J
    random(i,s)=randu(seed)
end do
end do

i=1
do s=1,typesN
    sim_types(1,i:i+int(simulationsN/(typesN))-1)=s
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(typesN))-1
end if
    i=i+int(simulationsN/(typesN))
end do


do s=2,J
    sim_types(s,:)=sim_types(1,:)
end do


seed=52345691
do s=1,simulationsN
do i=1,J
    random(i,s)=randu(seed)
end do
end do


i=1
do it = 1,typesN
    sim_shocks(1,i:i+int(simulationsN/(2*(typesN)))-1)=shocks_start1
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(2*(typesN)))-1
end if
    i=i+int(simulationsN/(2*(typesN)))
    sim_shocks(1,i:i+int(simulationsN/(2*(typesN)))-1)=shocks_start2
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(2*(typesN)))-1
end if
    i=i+int(simulationsN/(2*(typesN)))
end do

do i=1,shocksN
do s=1,shocksN
do it1=1,Jnr-1
do it2=1,2
    shock_cumpi(i,s,it1,it2)=sum(shock_pi(i,1:s,it1,it2))
end do
end do
end do
end do

do it1=Jnr,J
do it2=1,2
    shock_cumpi(:,:,it1,it2)=shock_cumpi(:,:,Jnr-1,it2)
end do
end do

shock_cumpi(:,shocksN,:,:)=1.0D0


!Low
do s=1,int(simulationsN/2)
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,1),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,1).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do
!High
do s=int(simulationsN/2)+1,simulationsN
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s)=minloc(shock_cumpi(sim_shocks(i-1,s),:,i,2),1,mask=shock_cumpi(sim_shocks(i-1,s),:,i,2).ge.random(i,s))
    else
        sim_shocks(i,s)=sim_shocks(i-1,s)
    end if
end do
end do

if(my_id .eq. root_process) then
print*,"I got herea"
end if
ss_count=0.0D0
do it =1 ,simulationsN
    if(sim_shocks(Jnr,it).le. shocksN/2) then
        ss_count(sim_types(1,it),sim_shocks(Jnr,it))=ss_count(sim_types(1,it),sim_shocks(Jnr,it))+1.0D0
    else
        ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)=ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)+1.0D0
    end if
end do


!MASTER
if (my_id .eq. root_process) then

do rebate_loop=1,slopeN*levelN*progressivityN+(1+progressivityN)*extra
    print*,"rebate loop",rebate_loop
if(routine==5) then
    rebate_level=rebate_level_ig
else
    if(vary_capital==1) then
        rebate_level=0.0D0
        lambdak=parameters(rebate_loop,1)
        rebate_slope=parameters(rebate_loop,2)
    elseif(vary_capital==2) then
        rebate_slope=0.0D0
        lambdak=parameters(rebate_loop,2)
        rebate_level=parameters(rebate_loop,1)
    else
        rebate_level=parameters(rebate_loop,1)
        rebate_slope=parameters(rebate_loop,2)
    end if
end if
lambda1=parameters(rebate_loop,3)+lambda1_benchmark
if(my_id .eq. root_process) then
print*,"I got here"
end if

if (routine ==4 .and. update_lambda0_ig==1 .and. rebate_loop>1) then
!offset other taxes
CALL lambda0_optimizer(lambda0_1,gov_spend-Ep*taue,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,&
types,partprod,r,tr,population(:,2),avg_earnings,ss,survive)
lambda0=lambda0_1
end if



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ENTER FUNCTION HERE IF MAKING A FUNCTION
!Initial guesses
K=K_ig
N=N_ig
tr=tr_ig
Ec=Ec_ig
ss=ss_ig
lump_sum=lump_sum_ig
if(my_id .eq. root_process) then
print*,"here3",rebate_level,rebate_slope,lambda1
end if

!!!!!!!!!!!!!!!!!!!!!!!! Initialize loop stuff
tax_loop = 1
tax_diff=1
non_clear=0
bad_combination_inner_stop=0
bad_combination_inner_nostop=0
bad_combination_outer=0
diff_parameter=1.0D0
out_loop=1
allocate(zeros2(max_outer_loops+1))
zeros2=0
taxes_clear=zeros2
deallocate(zeros2)
bad_combination_inner_stop=0
bad_combination_outer=0

!Loop for price convergence
DO WHILE (diff_parameter>tol)
   tax_loop = 1
   tax_diff=1.0D0
   non_clear=0

    passed_parameters(1)=K
    passed_parameters(2)=N
    passed_parameters(3)=Ec
    passed_parameters(4)=Atfp
    passed_parameters(5)=Ae
    passed_parameters(6)=taue
print*,"parameters_passed",passed_parameters

        allocate(share_guesses(6))
        allocate(share_guesses2(8))
        allocate(share_guesses3(8))
        if(out_loop==1 .and. rebate_loop==1) then
            share_guesses=(/.025D0,.035D0,.05D0,.015D0,.04D0,.01D0/)
        else
            share_guesses=(/.5D0,Ep,Ep*1.1D0,Ep*.8D0,Ep*1.05D0,Ep*.9D0/)
        end if
        share_guesses2=(/.5D0,.65D0,.8D0,.9D0,.94D0,.95D0,.96D0,.98D0/)
        share_guesses3=(/.5D0,.65D0,.8D0,.9D0,.94D0,.95D0,.96D0,.98D0/)
        do share_it=1,6
        do share_it2=1,8
        do share_it3=1,8
            xguesses_check(1)=share_guesses(share_it)
            xguesses_check(2)=K*share_guesses2(share_it2)
            xguesses_check(3)=N*share_guesses3(share_it3)
            CALL energy_root_finder(xguesses_check,check_dummy_out,3)
            if(share_it==1 .and. share_it2==1 .and. share_it3==1) then
                check_dummy=check_dummy_out(1)**2.0D0+check_dummy_out(2)**2.0D0+check_dummy_out(3)**2.0D0
                xguess(1)=xguesses_check(1)
                xguess(2)=xguesses_check(2)
                xguess(3)=xguesses_check(3)
            elseif (check_dummy>(check_dummy_out(1)**2.0D0+check_dummy_out(2)**2.0D0+check_dummy_out(3)**2.0D0)) then
                xguess(1)=xguesses_check(1)
                xguess(2)=xguesses_check(2)
                xguess(3)=xguesses_check(3)
                check_dummy=check_dummy_out(1)**2.0D0+check_dummy_out(2)**2.0D0+check_dummy_out(3)**2.0D0
            end if
        end do
        end do
        end do
        deallocate(share_guesses)
        deallocate(share_guesses2)
        deallocate(share_guesses3)
print*,"xguess in",xguess

    CALL D_NEQNF(energy_root_finder,solver_out,ERRREL=0.0000001D0,xguess=xguess)
print*,"xguess out",solver_out

    Ep=solver_out(1)
    E=Ep+Ec
    K1=solver_out(2)
    K2=K-K1
    N1=solver_out(3)
    N2=N-N1


    Y=Atfp*(K1**(alpha1)*N1**(1.0D0-alpha1-theta))*Ep**theta
    pe=Atfp*theta*(K1/Ep)**alpha1*(N1/Ep)**(1.0D0-alpha1-theta)-taue
    w=Atfp*(1.0D0-alpha1-theta)*(K1/N1)**(alpha1)*(Ep/N1)**theta
    r=Atfp*alpha1*(N1/K1)**(1.0D0-alpha1-theta)*(Ep/K1)**theta-delta
    tauss=sum(sum(ss*ss_count,2),1)/real(simulationsN)*sum(population(Jnr:J,2))/&
        (sum(sum(avg_earnings_ss*ss_count,2),1)/real(simulationsN)*sum(population(1:Jnr-1,2)))

!    tauss=sum(sum(ss*ss_count))/real(simulationsN)*sum(population(Jnr:J,2))/(avg_earnings_ss)
    ss_cap=avg_earnings*2.42
print*,"prices",w,r,avg_earnings_scale,pe,tauss,lambda0,lambdak

!Extra loop for convergence on tax parameters (if turned on need tax parameter convergence begore prices update)
    DO WHILE (tax_diff >tol)
!Solve for policy funcation
        v = -1.0D0/0.0D0
        v(:,:,:,J+1)=0.0D0
        v_retired = -1.0D0/0.0D0
        kprime=0
        v_retired(:,:,capital_positive:kgridN,J+1)=0.0D0

!Last generation
            do it=capital_positive,kgridN
                do it2=1,typesN
                    do e_it=1,shocksN
                        if(e_it .le. (shocksN)/2) then
                            ctot80(it,1)=retcon(kgrid(it),0.0D0,tr,r,ss(it2,e_it),avg_earnings)
                            energy_miss=energy_sharer(0.00D0,ctot80(it,1)*.15D0,ctot80(it,1),energy_share_equation,.0001D0,c80(it,1),ctot80(it,1),J)
                            e80(it,1)=(ctot80(it,1)-c80(it,1))/(pe+taue)
                            v_retired(it2,e_it,it,J) =utility_last(c80(it,1),e80(it,1),adult_equivalence(J))


                        else
                            v_retired(it2,e_it,it,J) =v_retired(it2,e_it-(shocksN/2),it,J)
                        end if
                    end do
                end do
            end do


!Retired generations
        allocate(zeros(kgridN,1))
        zeros=0.0D0
        do s =(J-1),Jnr, -1
            do e_it=1,shocksN
            if(e_it .le. (shocksN)/2) then
                do it=1,typesN
                do kit=capital_positive, kgridN
                    con=-1.0D0
                    con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,kgrid(kit),kgrid(capital_positive:),&
                        tr,r,ss(it,e_it),avg_earnings)
                    con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                    if (con(kgridN)>0.015+ebar) then
                        con_pos=kgridN
                    end if
                    if (con(capital_positive+1)<=0.015+ebar) then
                        kprime(it,e_it,kit,s)=0
                        v_retired(it,e_it,kit,s)=-1.0D0/0.0D0
                    else
                    CALL retired_valuer(kgrid(kit),0.0D0,kgrid(con_pos),ss(it,e_it),tr,r,survive(s),v_retired(it,e_it,:,s+1),&
                        v_retired(it,e_it,kit,s),kprime(it,e_it,kit,s),avg_earnings,s)
                    end if
               end do
               end do
            else
                do it=1,typesN
                do kit=capital_positive, kgridN
                    kprime(it,e_it,kit,s)=kprime(it,e_it-(shocksN/2),kit,s)
                    v_retired(it,e_it,kit,s)=v_retired(it,e_it-(shocksN/2),kit,s)
                end do
                end do
            end if
            end do
        end do
        deallocate(zeros)

if(my_id==1) then
print*,"here4"
end if


!Working generations
v(:,:,:,Jnr)=v_retired(:,:,:,Jnr)

                prices(1)=tr
                prices(2)=lambda0
                prices(3)=w
                prices(4)=r
                prices(5)=tauss
                if(routine==5) then
                    prices(6)=rebate_level
                end if
!                prices(6)=ss
                prices(7)=ss_cap
                prices(8)=lump_sum
                prices(9)=avg_earnings
                prices(10)=lambdak
                prices(11)=pe
                CALL MPI_Bcast( prices, 11, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
                CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
!print*,"prices sends",prices
        do s =(Jnr-1),2,-1
            CALL MPI_Bcast( v(:,:,:,s+1),typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            end_row=0
            do an_id = 1, num_procs_run -1
                start_row=end_row+1
                if (an_id>large_inc) then
                    end_row=start_row+increment-1
                else
                    end_row=start_row+increment
                end if
                call MPI_SEND(start_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
                call MPI_SEND(end_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
            end do
            do an_id=1,num_procs_run -1
                call MPI_RECV( start_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( end_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( labor(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( kprime(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( v(:,:,start_row:end_row,s),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            end do

        end do

if(my_id==1) then
print*,"here5"
end if




!! Simulate peoples

    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end_row=0
    do an_id = 1, num_procs_run2 -1
        start_row=end_row+1
        if (an_id>large_inc2) then
            end_row=start_row+increment2-1
        else
            end_row=start_row+increment2
        end if
        call MPI_SEND(start_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
        call MPI_SEND(end_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
    end do

    do an_id=1,num_procs_run2 -1
        call MPI_RECV( start_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( end_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_labor(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_earnings(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_capital(:,start_row:end_row),(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_consumption(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_energy(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_v(start_row:end_row),(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_total_income(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_after_tax_income(:,start_row:end_row),(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    end do

print*,"value",sum(sim_v)/real(simulationsN)



!Calculate Aggregate Profiles

agg_labor(:,1)=sum(sim_labor(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_labor(:,2)=sum(sim_labor(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_capital(:,1)=sum(sim_capital(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_capital(:,2)=sum(sim_capital(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_consumption(:,1)=sum(sim_consumption(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_consumption(:,2)=sum(sim_consumption(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_energy(:,1)=sum(sim_energy(:,1:int((simulationsN)/2)),2)/int((simulationsN)/2.0D0)
agg_energy(:,2)=sum(sim_energy(:,int((simulationsN)/2)+1:simulationsN),2)/(simulationsN-int((simulationsN)/2.0D0))





ind_eng_tax=0.0D0
do an_id=1,J
    ind_eng_tax=ind_eng_tax+taue*(agg_energy(an_id,1)*.5D0+.5D0*agg_energy(an_id,2))*population(an_id,2)
end do
!!Find lambda 0 that clears the market


!1 lump_sum_rebate
!2 throw_away
!3 cap_tax_rebate
!4 labor_rebate

if (routine == 3) then
CALL lambdak_optimizer(lambda0_1,gov_spend-Ep*taue,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,&
                    types,partprod,r,tr,population(:,2),avg_earnings,ss,survive)




elseif (routine ==4) then
!offset other taxes
CALL lambda0_optimizer(lambda0_1,gov_spend-Ep*taue,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,&
                    types,partprod,r,tr,population(:,2),avg_earnings,ss,survive)

elseif (routine ==5) then
!offset other taxes
CALL rebate_optimizer(rebate_level1,gov_spend-Ep*taue,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,&
                    types,partprod,r,tr,population(:,2),avg_earnings,ss,survive)!rebate



else if( routine==1  ) then
allocate(dummy_zeros(J,simulationsN))
dummy_zeros=0.0D0
CALL lambda0_optimizer(lambda0_1,gov_spend,sim_labor,dummy_zeros,sim_shocks,sim_types,sim_capital,hbar,w,shocks,&
                    types,partprod,r,tr,population(:,2),avg_earnings,ss,survive)
else if(routine==2) then
    lambda0_1=lambda0
else
    print*,"Not Specified Right"
end if


        if (routine == 3) then
            tax_diff=abs((lambda0_1)-lambdak)/lambdak
            tax_diff_keep=tax_diff
print*,"tax diff",tax_diff,lambda0_1,lambdak
            lambdak=(a1*lambda0_1+(1.0D0-a1)*lambdak)
        elseif(routine==5) then
            tax_diff=abs((rebate_level1)-rebate_level)/rebate_level
            tax_diff_keep=tax_diff
print*,"tax diff",tax_diff,rebate_level1,rebate_level
            rebate_level=(a1*rebate_level1+(1.0D0-a1)*rebate_level)
        elseif(routine==4) then
        govt_taxes1=tax_differ2(lambda0,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,gov_spend-Ep*taue,population(:,2),avg_earnings,ss,lambdak,survive,rebate_level)
        govt_taxes2=tax_differ2(lambda0_1,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,gov_spend-Ep*taue,population(:,2),avg_earnings,ss,lambdak,survive,rebate_level)
        govt_taxes3=tax_differ2(0.4D0,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,gov_spend-Ep*taue,population(:,2),avg_earnings,ss,lambdak,survive,rebate_level)
        govt_taxes4=tax_differ2(.90D0,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,gov_spend-Ep*taue,population(:,2),avg_earnings,ss,lambdak,survive,rebate_level)
        if(abs(govt_taxes1)<abs(abs(govt_taxes2)))then
            lambda0_1=lambda0
        end if

            tax_diff=abs((lambda0_1)-lambda0)/lambda0
            tax_diff_keep=tax_diff
print*,"tax diff",tax_diff,lambda0_1,lambda0,gov_spend-Ep*taue,govt_taxes1,govt_taxes2,govt_taxes3,govt_taxes4
        lambda0=(a1*lambda0_1+(1.0D0-a1)*lambda0)
        elseif(routine==2) then
            tax_diff=0.0D0
            tax_diff_keep=tax_diff
        else
            tax_diff=abs((lambda0_1)-lambda0)/lambda0
            tax_diff_keep=tax_diff
            print*,"tax diff",tax_diff,lambda0_1,lambda0,gov_spend-Ep*taue
            lambda0=(a1*lambda0_1+(1.0D0-a1)*lambda0)
        end if

if (routine == 3 .or. routine ==4 .or. routine==5) then
!offset other taxes
        govt_taxes=tax_differ2(lambda0,sim_labor,sim_energy,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,gov_spend-Ep*taue,population(:,2),avg_earnings,ss,lambdak,survive,rebate_level)
else if( routine==1  ) then
!rebate
        govt_taxes=tax_differ2(lambda0,sim_labor,dummy_zeros,sim_shocks,sim_types,sim_capital,hbar,w,shocks,types,partprod,&
                    r,tr,gov_spend,population(:,2),avg_earnings,ss,lambdak,survive,rebate_level)
deallocate(dummy_zeros)
elseif(routine==2) then
    govt_taxes=0.0D0

else
    print*,"Not Specified Right"
end if

        !Check for bad convergence within tax loop

        tax_loop=tax_loop+1
        if (tax_loop > max_tax_loop) then
          tax_diff = 0
          bad_combination_inner_nostop=1
        end if
        CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    END DO





    K1_prime=0.0D0
    N1_prime=0.0D0
    tr1=0.0D0
    hours_jnr=0.0D0
    Ec1=0.0D0

    avg_earnings1=0.0D0
    avg_earnings_ss1=0.0D0


    !Calculate new aggregates
    do it=1,simulationsN
        do s=1,Jnr-1
            avg_earnings1=avg_earnings1+sim_earnings(s,it)*population(s,2)/(real(simulationsN)*sum(population(1:Jnr-1,2)))
            if (sim_shocks(Jnr,it) .le. shocksN/2) then
                avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it))=&
                    avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it))+&
                    min(ss_cap,sim_earnings(s,it))*population(s,2)/(ss_count(sim_types(1,it),sim_shocks(Jnr,it))*&
                    sum(population(1:Jnr-1,2)))
            else
                avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)=&
                    avg_earnings_ss1(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)+&
                    min(ss_cap,sim_earnings(s,it))*population(s,2)/(ss_count(sim_types(1,it),sim_shocks(Jnr,it)-shocksN/2)*&
                    sum(population(1:Jnr-1,2)))
            end if
            Ec1=Ec1+sim_energy(s,it)*population(s,2)/(simulationsN)
        end do
        do s=1,Jnr-1
            if (sim_labor(s,it)<hbar) then
                N1_prime=N1_prime+partprod*sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN
            else
                N1_prime=N1_prime+sim_labor(s,it)*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*population(s,2)/simulationsN
            end if
            K1_prime=K1_prime+sim_capital(s+1,it)*population(s,2)/simulationsN
            tr1=tr1+sim_capital(s+1,it)*population(s,2)/simulationsN*(1.0D0-survive(s))
            hours_jnr=hours_jnr+sim_labor(s,it)*population(s,2)/simulationsN
        end do
        do s=Jnr,J
            Ec1=Ec1+sim_energy(s,it)*population(s,2)/(simulationsN)
            K1_prime=K1_prime+sim_capital(s+1,it)*population(s,2)/simulationsN
            tr1=tr1+sim_capital(s+1,it)*population(s,2)/simulationsN*(1.0D0-survive(s))
        end do
    end do


!Lump sum rebate
if(routine==1) then
    lump_sum=ind_eng_tax+taue*E
elseif (routine==2 .or. routine ==3 .or. routine ==4) then
    lump_sum=0.0D0
else
    print*,"Not Specified Right"
end if

        do e_it=1,shocksN/2
            do it=1,typesN
                ss1(it,e_it)=ss_start(it,e_it)*(Ce_for_share_start*(taue+pe)+C_for_share_start)/(Ce_for_share_start*pe_start+C_for_share_start)
            end do
        end do




        hours_jnr=hours_jnr/sum(population(1:Jnr-1,2))
        K1_prime=K1_prime/(1.0D0+growth)
        tr1=tr1/(1.0D0+growth)

    diff_parameter_v(1)=abs(((tr1)-tr))/tr
    diff_parameter_v(2)=abs(((K1_prime)-K))/K
    diff_parameter_v(3)=abs(((N1_prime)-N))/N
!    diff_parameter_v(4)=abs(((ss1*a+ss*(1.0D0-a))-ss))/ss
    diff_parameter_v(4)=maxval(abs(((ss1)-ss))/ss)
!    ,abs(((ss1(2))-ss(2)))/ss(2))
    diff_parameter_v(5)=abs(tax_diff_keep)
    diff_parameter_v(6)=max(abs(((avg_earnings1)-avg_earnings))/avg_earnings,abs(((Ec1)-Ec))/Ec&
        ,maxval(abs(((avg_earnings_ss1)-avg_earnings_ss))/avg_earnings_ss1))
    diff_parameter=MAXVAL(diff_parameter_v)
    print*,"K1",K1_prime,K
    print*,"N1",N1_prime,N
    print*,"N1",Ec1,Ec
    print*,"tr1",tr1,tr
    print*,"ss1",ss1,ss
    print*,"avg_earnings1",avg_earnings1,avg_earnings
!    print*,"avg_earnings_ss1",avg_earnings_ss1,avg_earnings_ss
    print*,"hours",hours_jnr
    print*,"lump sum",lump_sum

!!!!Update Aggregates
    Ec=(Ec1*a+Ec*(1-a))
    N=(N1_prime*a+N*(1-a))
    K=(K1_prime*a+K*(1-a))
    ss=(ss1*a+ss*(1-a))
    tr=(tr1*a+tr*(1-a))
    avg_earnings=(avg_earnings1*a+avg_earnings*(1-a))
    avg_earnings_ss=(avg_earnings_ss1*a+avg_earnings_ss*(1-a))

    if (bad_combination_inner_stop == 1) then
        diff_parameter=0.0D0
    end if
    taxes_clear(out_loop)=0
    if (abs(govt_taxes/gov_spend)>tol) then
        taxes_clear(out_loop)=1
        non_clear=1
        diff_parameter=1.0D0
    end if
    if (taxes_clear(out_loop)==0) then
        allocate(zeros2(max_outer_loops+1))
        zeros2=0
        taxes_clear=zeros2
        deallocate(zeros2)
    end if
    if (sum(taxes_clear)>max_non_clear) then
        diff_parameter=0.0D0
        bad_combination_inner_stop=1
    end if
    if (out_loop > max_outer_loops) then
        bad_combination_outer=1
        diff_parameter=0.0D0
    end if
!end if
    out_loop=out_loop+1
    print*,"Loop"
    print*,"Plus one"
    print*,rebate_loop,out_loop
    print*,"lambda0",lambda0
    print*,"lambdak",lambdak
    print*,"govt miss",govt_taxes/gov_spend
    print*,"diffparm",diff_parameter
    print*,"taxdiff",tax_diff_keep
    print*,"hours",hours_jnr
    print*,"K",K
    print*,"N",N
    print*,"Y",Y
    print*,"ss",ss(1,:)
    print*,"ss",ss(2,:)
    print*,"tr",tr
    print*,"E",E
    print*,"w",w
    print*,"r",r
    print*,"income", w*shocks(sim_shocks(2,1))*types(sim_types(2,1))*epsilon1(2,1),sim_shocks(2,1)
    print*,"income2", w*shocks(sim_shocks(3,1))*types(sim_types(3,1))*epsilon1(3,1),sim_shocks(3,1)
    print*,"choices",sim_earnings(1,1),sim_capital(1,1),sim_consumption(1,1),sim_energy(1,1),sim_total_income(1,1),sim_after_tax_income(1,1)
    print*,"choices",sim_earnings(2,1),sim_capital(2,1),sim_consumption(2,1),sim_energy(2,1),sim_total_income(2,1),sim_after_tax_income(2,1)
    print*,"choices",sim_earnings(3,1),sim_capital(3,1),sim_consumption(3,1),sim_energy(3,1),sim_total_income(3,1),sim_after_tax_income(3,1)
    print*,"choices",sim_earnings(4,1),sim_capital(4,1),sim_consumption(4,1),sim_energy(4,1),sim_total_income(4,1),sim_after_tax_income(4,1)

    print*,bad_combination_inner_nostop
    CALL MPI_Bcast( diff_parameter, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )

END DO

sim_lifetime_total_income=0.0D0
sim_lifetime_after_tax_income=0.0D0
sim_lifetime_total_income_nosurvive=0.0D0
sim_lifetime_after_tax_income_nosurvive=0.0D0

do it =J,1,-1
    do it2=1,simulationsN
        ypp_temp=sim_earnings(it,it2)-.5D0*tauss*min(sim_earnings(it,it2),ss_cap)
        sim_labor_tax_new(it,it2)=lab_taxer(ypp_temp)
        sim_labor_tax_old(it,it2)=old_lab_taxer(ypp_temp)
        if(it<Jnr-1) then
            sim_rebate(it,it2)=rebate_calculator(sim_earnings(it,it2),(sim_capital(it,it2)+tr)*(r),&
            1,sim_labor_tax_new(it,it2),lambdak*(sim_capital(it,it2)+tr)*(r))
        else
            sim_rebate(it,it2)=rebate_calculator(0.0D0,(sim_capital(it,it2)+tr)*(r),&
            0,0.0D0,lambdak*(sim_capital(it,it2)+tr)*(r))
        end if
    end do
end do


do it =J,1,-1
    do it2=1,simulationsN
        if (it==J) then
            sim_utility(it,it2)=utility_last(sim_consumption(it,it2),sim_energy(it,it2),adult_equivalence(it))
            sim_utility_start(it,it2)=utility_last(sim_consumption_start(it,it2),sim_energy_start(it,it2),adult_equivalence(it))

            sim_lifetime_total_income(it2)=sim_total_income(it,it2)
            sim_lifetime_after_tax_income(it2)=sim_after_tax_income(it,it2)
            sim_lifetime_after_tax_income_start(it2)=sim_total_after_tax_earnings_start(it,it2)


            sim_lifetime_total_income_nosurvive(it2)=sim_total_income(it,it2)
            sim_lifetime_after_tax_income_nosurvive(it2)=sim_after_tax_income(it,it2)
        elseif(it>Jnr-1) then
            sim_utility(it,it2)=valuer_ret(sim_consumption(it,it2),sim_energy(it,it2),survive(it),sim_utility(it+1,it2),adult_equivalence(it))
            sim_utility_start(it,it2)=valuer_ret(sim_consumption_start(it,it2),sim_energy_start(it,it2),survive(it),sim_utility_start(it+1,it2),adult_equivalence(it))


            sim_lifetime_total_income(it2)=survive(it)*1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_total_income(it2)+sim_total_income(it,it2)
            sim_lifetime_after_tax_income(it2)=survive(it)*1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_after_tax_income(it2)+sim_after_tax_income(it,it2)
            sim_lifetime_after_tax_income_start(it2)=1.0D0/(1.0D0+r_start*(1-lambdak_start))*sim_lifetime_after_tax_income_start(it)+sim_total_after_tax_earnings_start(it,it2)

            sim_lifetime_total_income_nosurvive(it2)=1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_total_income_nosurvive(it2)+sim_total_income(it,it2)
            sim_lifetime_after_tax_income_nosurvive(it2)=1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_after_tax_income_nosurvive(it2)+sim_after_tax_income(it,it2)
        else
            sim_utility(it,it2)=valuer(sim_consumption(it,it2),sim_energy(it,it2),sim_labor(it,it2),survive(it),sim_utility(it+1,it2),adult_equivalence(it))
            sim_utility_start(it,it2)=valuer(sim_consumption_start(it,it2),sim_energy_start(it,it2),sim_labor_start(it,it2),survive(it),sim_utility_start(it+1,it2),adult_equivalence(it))

            sim_lifetime_total_income(it2)=survive(it)*1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_total_income(it2)+sim_total_income(it,it2)
            sim_lifetime_after_tax_income(it2)=survive(it)*1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_after_tax_income(it2)+sim_after_tax_income(it,it2)
            sim_lifetime_after_tax_income_start(it2)=1.0D0/(1.0D0+r_start*(1-lambdak_start))*sim_lifetime_after_tax_income_start(it)+sim_total_after_tax_earnings_start(it,it2)

            sim_lifetime_total_income_nosurvive(it2)=1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_total_income_nosurvive(it2)+sim_total_income(it,it2)
            sim_lifetime_after_tax_income_nosurvive(it2)=1.0D0/(1.0D0+r*(1-lambdak))*sim_lifetime_after_tax_income_nosurvive(it2)+sim_after_tax_income(it,it2)
        end if
        sim_total_consumption_start(it,it2)=sim_consumption_start(it,it2)+sim_energy_start(it,it2)*pe
        sim_consumption_tilda(it,it2)=consumption_tilda(sim_consumption(it,it2),sim_energy(it,it2))
    end do
end do


do it =J,1,-1
    do it2=1,simulationsN
        if (it==J) then
            sim_consumption_dummy(1,it2)=sim_consumption(it,it2)
            sim_consumption_dummy_start(1,it2)=sim_consumption_start(it,it2)


        elseif(it>Jnr-1) then
            sim_consumption_dummy(1,it2)=sim_consumption(it,it2)+survive(it)*beta*sim_consumption_dummy(1,it2)
            sim_consumption_dummy_start(1,it2)=sim_consumption_start(it,it2)+survive(it)*beta*sim_consumption_dummy_start(1,it2)


        else
            sim_consumption_dummy(1,it2)=sim_consumption(it,it2)+survive(it)*beta*sim_consumption_dummy(1,it2)
            sim_consumption_dummy_start(1,it2)=sim_consumption_start(it,it2)+survive(it)*beta*sim_consumption_dummy_start(1,it2)

            sim_labor_dummy(1,it2)=sim_labor(it,it2)+survive(it)*beta*sim_labor_dummy(1,it2)
            sim_labor_dummy_start(1,it2)=sim_labor_start(it,it2)+survive(it)*beta*sim_labor_dummy_start(1,it2)

        end if
    end do
end do



perc_better=0.0D0
do it=1,simulationsN
    if(sim_utility(1,it)>sim_utility_start(1,it)) then
        perc_better=perc_better+1.0D0
    else
    end if
end do
perc_better=perc_better/real(simulationsN)
!print*,"better off",perc_better

!Gini for Welfare
giniWelSS_start=0.0D0
lifeWelBar=sum(sim_utility_start(1,:))/real(simulationsN)
giniadj=abs(2.0D0*real(simulationsN)**2.0D0*lifeWelBar)
do it =1,simulationsN
    do it2=1,simulationsN
        giniWelSS_start=giniWelSS_start+abs(sim_utility_start(1,it)-sim_utility_start(1,it2))/giniadj
    end do
end do

!Gini for Welfare
giniWelSS=0.0D0
lifeWelBar=sum(sim_utility(1,:))/real(simulationsN)
giniadj=abs(2.0D0*real(simulationsN)**2.0D0*lifeWelBar)
do it =1,simulationsN
    do it2=1,simulationsN
        giniWelSS=giniWelSS+abs(sim_utility(1,it)-sim_utility(1,it2))/giniadj
    end do
end do
!print*,"gini welfare",giniWelSS,giniWelSS_start


!Gini for after tax income
giniIncSS_aftertax=0.0D0
lifeWelBar=sum(sim_lifetime_after_tax_income(:))/real(simulationsN)
giniadj=abs(2.0D0*real(simulationsN)**2.0D0*lifeWelBar)
do it =1,simulationsN
    do it2=1,simulationsN
        giniIncSS_aftertax=giniIncSS_aftertax+abs(sim_lifetime_after_tax_income(it)-sim_lifetime_after_tax_income(it2))/giniadj
    end do
end do

!Gini for total tax income
giniIncSS=0.0D0
lifeWelBar=sum(sim_lifetime_total_income(:))/real(simulationsN)
giniadj=abs(2.0D0*real(simulationsN)**2.0D0*lifeWelBar)
do it =1,simulationsN
    do it2=1,simulationsN
        giniIncSS=giniIncSS+abs(sim_lifetime_total_income(it)-sim_lifetime_total_income(it2))/giniadj
    end do
end do


!Gini for after tax income
giniIncSS_aftertax_nosurvive=0.0D0
lifeWelBar=sum(sim_lifetime_after_tax_income_nosurvive(:))/real(simulationsN)
giniadj=abs(2.0D0*real(simulationsN)**2.0D0*lifeWelBar)
do it =1,simulationsN
    do it2=1,simulationsN
        giniIncSS_aftertax_nosurvive=giniIncSS_aftertax_nosurvive+abs(sim_lifetime_after_tax_income_nosurvive(it)-sim_lifetime_after_tax_income_nosurvive(it2))/giniadj
    end do
end do

!Gini for total tax income
giniIncSS_nosurvive=0.0D0
lifeWelBar=sum(sim_lifetime_total_income_nosurvive(:))/real(simulationsN)
giniadj=abs(2.0D0*real(simulationsN)**2.0D0*lifeWelBar)
do it =1,simulationsN
    do it2=1,simulationsN
        giniIncSS_nosurvive=giniIncSS_nosurvive+abs(sim_lifetime_total_income_nosurvive(it)-sim_lifetime_total_income_nosurvive(it2))/giniadj
    end do
end do





CALL cev_finder(cev_total,sim_consumption_start,sim_energy_start,sim_labor_start,survive,&
    sim_consumption,sim_energy,sim_labor,simulationsN,0)

CALL cev_finder(cev_total_alt,sim_consumption_start,sim_energy_start,sim_labor_start,survive,&
    sim_consumption,sim_energy,sim_labor,simulationsN,1)


CALL cev_finder(cev_total_collapsed,sum(sim_consumption_start,2)/real(simulationsN),sum(sim_energy_start,2)/real(simulationsN),sum(sim_labor_start,2)/real(simulationsN),survive,&
    sum(sim_consumption,2)/real(simulationsN),sum(sim_energy,2)/real(simulationsN),sum(sim_labor,2)/real(simulationsN),1,0)



!Output for steady state to be used in transition (saved before sorting)
!         OPEN(UNIT=14, FILE="/<filepath steady state for transition output>/v_retired_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!          WRITE(14,887) v_retired
!          CLOSE(UNIT=14)
!
!
!         OPEN(UNIT=15, FILE="/<filepath steady state for transition output>/sim_labor_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!          do it=1,simulationsN
!          do it2=1,J
!          WRITE(15,887) sim_labor(it2,it)
!            end do
!            end do
!          CLOSE(UNIT=15)
!
!          OPEN(UNIT=16, FILE="/<filepath steady state for transition output>/sim_energy_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!          do it=1,simulationsN
!          do it2=1,J
!          WRITE(16,887) sim_energy(it2,it)
!            end do
!            end do
!          CLOSE(UNIT=16)
!
!         OPEN(UNIT=17, FILE="/<filepath steady state for transition output>/sim_capital_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!          do it=1,simulationsN
!          do it2=1,J+1
!          WRITE(17,887) sim_capital(it2,it)
!            end do
!            end do
!          CLOSE(UNIT=17)
!
!         OPEN(UNIT=18, FILE="/<filepath steady state for transition output>/sim_consumption_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!          do it=1,simulationsN
!          do it2=1,J
!          WRITE(18,887) sim_consumption(it2,it)
!            end do
!            end do
!          CLOSE(UNIT=18)
!
!         OPEN(UNIT=13, FILE="/<filepath steady state for transition output>/v_output.dat",STATUS="replace", &
!               FORM="formatted")
!          WRITE(13,887) v
!          CLOSE(UNIT=13)
!!
!!         OPEN(UNIT=19, FILE="/<filepath steady state for transition output>/optimal.dat",STATUS="replace", &
!!               FORM="formatted")
!!          WRITE(19,887) aggregates(:,rebate_loop)
!!          CLOSE(UNIT=13)
!
!        OPEN(UNIT=21, FILE="/<filepath steady state for transition output>/sim_shocks_output.dat",STATUS="replace", &
!               FORM="formatted")
!          do it=1,simulationsN
!          do it2=1,J
!          WRITE(21,886) sim_shocks(it2,it)
!            end do
!            end do
!          CLOSE(UNIT=21)
!
!         OPEN(UNIT=15, FILE="/<filepath steady state for transition output>/ss_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!          WRITE(15,887) ss
!          CLOSE(UNIT=15)
!
!         OPEN(UNIT=15, FILE="/<filepath steady state for transition output>/avg_earnings_ss_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!          WRITE(15,887) avg_earnings_ss
!          CLOSE(UNIT=15)
!
!
!         OPEN(UNIT=15, FILE="/<filepath steady state for transition output>/avg_earnings_ss_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!          WRITE(15,887) avg_earnings_ss
!          CLOSE(UNIT=15)
!
!         OPEN(UNIT=18, FILE="/<transition output path>/sim_rebate_output_kevin.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!               do it=1,simulationsN
!          WRITE(18,887) sim_rebate(:,it)
!          end do
!          CLOSE(UNIT=18)




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Sort and CEV by lifetime utility
    do it= 1, simulationsN
        krow = minloc( sim_utility_start( 1,it:simulationsN ), dim=1 ) + it - 1
        buf( : )     = sim_utility_start( :,it )
        sim_utility_start( :,it ) = sim_utility_start( :,krow )
        sim_utility_start( :,krow ) = buf( : )

        buf( : )     = sim_utility( :,it )
        sim_utility( :,it ) = sim_utility( :,krow )
        sim_utility( :,krow ) = buf( : )

        buf( : )     = sim_consumption( :,it )
        sim_consumption( :,it ) = sim_consumption( :,krow )
        sim_consumption( :,krow ) = buf( : )

        buf( : )     = sim_energy( :,it )
        sim_energy( :,it ) = sim_energy( :,krow )
        sim_energy( :,krow ) = buf( : )

        buf( : )     = sim_labor( :,it )
        sim_labor( :,it ) = sim_labor( :,krow )
        sim_labor( :,krow ) = buf( : )

        buf( : )     = sim_labor_start( :,it )
        sim_labor_start( :,it ) = sim_labor_start( :,krow )
        sim_labor_start( :,krow ) = buf( : )

        buf( : )     = sim_energy_start( :,it )
        sim_energy_start( :,it ) = sim_energy_start( :,krow )
        sim_energy_start( :,krow ) = buf( : )

        buf( : )     = sim_consumption_start( :,it )
        sim_consumption_start( :,it ) = sim_consumption_start( :,krow )
        sim_consumption_start( :,krow ) = buf( : )

        buf( 1 )     = sim_lifetime_after_tax_income_start( it )
        sim_lifetime_after_tax_income_start( it ) = sim_lifetime_after_tax_income_start( krow )
        sim_lifetime_after_tax_income_start( krow ) = buf( 1 )
    end do







!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Sort and CEV by lifetime utility
    do it= 1, simulationsN
        krow = minloc( sim_lifetime_after_tax_income_start( it:simulationsN ), dim=1 ) + it - 1
        buf( 1 )     = sim_lifetime_after_tax_income_start( it )
        sim_lifetime_after_tax_income_start( it ) = sim_lifetime_after_tax_income_start( krow )
        sim_lifetime_after_tax_income_start( krow ) = buf( 1 )

        buf( : )     = sim_utility( :,it )
        sim_utility( :,it ) = sim_utility( :,krow )
        sim_utility( :,krow ) = buf( : )

        buf( : )     = sim_consumption( :,it )
        sim_consumption( :,it ) = sim_consumption( :,krow )
        sim_consumption( :,krow ) = buf( : )

        buf( : )     = sim_energy( :,it )
        sim_energy( :,it ) = sim_energy( :,krow )
        sim_energy( :,krow ) = buf( : )

        buf( : )     = sim_labor( :,it )
        sim_labor( :,it ) = sim_labor( :,krow )
        sim_labor( :,krow ) = buf( : )

        buf( : )     = sim_labor_start( :,it )
        sim_labor_start( :,it ) = sim_labor_start( :,krow )
        sim_labor_start( :,krow ) = buf( : )

        buf( : )     = sim_energy_start( :,it )
        sim_energy_start( :,it ) = sim_energy_start( :,krow )
        sim_energy_start( :,krow ) = buf( : )

        buf( : )     = sim_consumption_start( :,it )
        sim_consumption_start( :,it ) = sim_consumption_start( :,krow )
        sim_consumption_start( :,krow ) = buf( : )

        buf( 1 )     = sim_lifetime_after_tax_income_start( it )
        sim_lifetime_after_tax_income_start( it ) = sim_lifetime_after_tax_income_start( krow )
        sim_lifetime_after_tax_income_start( krow ) = buf( 1 )

        buf( : )     = sim_utility_start( :,it )
        sim_utility_start( :,it ) = sim_utility_start( :,krow )
        sim_utility_start( :,krow ) = buf( : )

    end do













!!!quintile1
C_end=0.0D0
Ce_end=0.0D0
C_start=0.0D0
Ce_start=0.0D0
do it=1,simulationsN
    do s=1,J
        Ce_end=Ce_end+sim_energy(s,it)*population(s,2)
    end do
end do

sim_lifetime_after_tax_income_sorted(1,:)=sim_lifetime_after_tax_income
do it= 1, simulationsN
    krow = minloc( sim_lifetime_after_tax_income_sorted( 1,it:simulationsN ), dim=1) + it - 1

    buf( 1 )     = sim_lifetime_after_tax_income_sorted( 1,it )
    sim_lifetime_after_tax_income_sorted( 1,it ) = sim_lifetime_after_tax_income_sorted( 1,krow )
    sim_lifetime_after_tax_income_sorted( 1,krow ) = buf( 1 )

end do

aggregates(:,rebate_loop)=0.0D00
aggregates(1,rebate_loop)=w
aggregates(2,rebate_loop)=r
aggregates(3,rebate_loop)=K
aggregates(4,rebate_loop)=N
aggregates(5,rebate_loop)=tr
aggregates(6,rebate_loop)=Ce_end
aggregates(7,rebate_loop)=giniWelSS
aggregates(8,rebate_loop)=gov_spend
aggregates(9,rebate_loop)=tauss
aggregates(10,rebate_loop)=no_carbon_tax
aggregates(11,rebate_loop)=E
aggregates(12,rebate_loop)=lambdak
aggregates(13,rebate_loop)=lambda0
aggregates(14,rebate_loop)=lambda1
aggregates(15,rebate_loop)=avg_earnings
aggregates(16,rebate_loop)=rebate_level
aggregates(17,rebate_loop)=rebate_slope
aggregates(18,rebate_loop)=sum(sim_v)/real(simulationsN)
aggregates(19,rebate_loop)=giniIncSS
aggregates(20,rebate_loop)=sim_lifetime_after_tax_income_sorted(1,1)
aggregates(21,rebate_loop)=sim_lifetime_after_tax_income_sorted(1,1+INT(simulationsN/5))
aggregates(22,rebate_loop)=sim_lifetime_after_tax_income_sorted(1,1+2*INT(simulationsN/5))
aggregates(23,rebate_loop)=sim_lifetime_after_tax_income_sorted(1,1+3*INT(simulationsN/5))
aggregates(24,rebate_loop)=sim_lifetime_after_tax_income_sorted(1,1+4*INT(simulationsN/5))
aggregates(25,rebate_loop)=sim_lifetime_after_tax_income_sorted(1,simulationsN)
aggregates(26,rebate_loop)=cev_total
aggregates(50,rebate_loop)=perc_better
aggregates(56,rebate_loop)=out_loop
aggregates(57,rebate_loop)=(Ce_end+E)*taue
aggregates(58,rebate_loop)=hours_jnr
aggregates(64,rebate_loop)=0
aggregates(79,rebate_loop)=Ec
aggregates(80,rebate_loop)=K1
aggregates(81,rebate_loop)=N1
aggregates(82,rebate_loop)=K2
aggregates(83,rebate_loop)=N2
aggregates(84,rebate_loop)=Ep
aggregates(85,rebate_loop)=Pe
aggregates(86,rebate_loop)=sum(sum(sim_consumption(:,:),2)*population(1:J,2))
aggregates(87,rebate_loop)=sum(sum(sim_consumption_tilda(:,:),2)*population(1:J,2))
aggregates(88,rebate_loop)=sum(sum(sim_energy(:,:),2)*population(1:J,2))
aggregates(89,rebate_loop)=giniIncSS_aftertax
aggregates(90,rebate_loop)=giniIncSS_aftertax_nosurvive
aggregates(91,rebate_loop)=giniIncSS_nosurvive
aggregates(92,rebate_loop)=govt_taxes/gov_spend
aggregates(93,rebate_loop)=maxval(sim_capital)

aggregates(94,rebate_loop)=routine
aggregates(95,rebate_loop)=labor_capital_rebate_base
aggregates(96,rebate_loop)=vary_capital
aggregates(97,rebate_loop)=no_increase
aggregates(98,rebate_loop)=labor_capital_rebate_base
aggregates(99,rebate_loop)=rebate_negative_tax
aggregates(100,rebate_loop)=negative_taxes
aggregates(101,rebate_loop)=rebate_to_old
aggregates(102,rebate_loop)=flat_tax
aggregates(103,rebate_loop)=bad_combination_outer
aggregates(108,rebate_loop)=cev_total_alt


print*,"out",aggregates(12:14,rebate_loop),aggregates(16:17,rebate_loop),aggregates(26:29,rebate_loop)
!if(out_loop<max_outer_loops) then
    Ec_ig=Ec
    K_ig=K
    N_ig=N
!    E_ig=E
    tr_ig=tr
    ss_ig=ss
    lump_sum_ig=lump_sum
if(routine==5) then
    rebate_level_ig=rebate_level
end if


!Rebate loop
end do

!1 lump_sum_rebate
!2 throw_away
!3 cap_tax_rebate
!4 labor_rebate
!if(routine==3) then

!Output for search of optimal
         OPEN(UNIT=19, FILE="/<filepath optimal search output>//unconstrained_optimal_baseline2.dat",STATUS="replace", &
               FORM="formatted")
      do it2=1,slopeN*levelN*progressivityN+(progressivityN+1)*extra
          WRITE(19,987) aggregates(:,it2)
          end do

!Output for profiles to examine a steady state (these are post sorted so will not work for the steady state starting equilibrium for transition)
!         OPEN(UNIT=15, FILE="/<filepath optimal search output>//sim_labor_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!               do it=1,simulationsN
!          WRITE(15,887) sim_labor(:,it)
!          end do
!          CLOSE(UNIT=15)
!
!          OPEN(UNIT=16, FILE="/<filepath optimal search output>//sim_energy_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!               do it=1,simulationsN
!          WRITE(16,887) sim_energy(:,it)
!          end do
!          CLOSE(UNIT=16)
!
!         OPEN(UNIT=17, FILE="/<filepath optimal search output>//sim_capital_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!               do it=1,simulationsN
!          WRITE(17,887) sim_capital(:,it)
!          end do
!          CLOSE(UNIT=17)
!
!         OPEN(UNIT=18, FILE="/<filepath optimal search output>//sim_consumption_output.dat", ACTION="write", STATUS="replace", &
!               FORM="formatted")
!               do it=1,simulationsN
!          WRITE(18,887) sim_consumption(:,it)
!          end do
!          CLOSE(UNIT=18)
!
!
!!


987 format(135F24.16)

886 format(i9)
887 format(F24.16)
889 format(2F24.16)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Slaves
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!Slaves that are going to work
else if (my_id < (num_procs_run)) then
do rebate_loop=1,slopeN*levelN*progressivityN+(1+progressivityN)*extra
!!!!!!!!!!!!!!!!!!!!!!!! Initialize loop stuff
out_loop=1
tax_loop = 1
tax_diff=1
non_clear=0
bad_combination_inner_stop=0
bad_combination_inner_nostop=0
bad_combination_outer=0
diff_parameter=1.0D0

if(routine==5) then
    rebate_level=rebate_level_ig
else
    if(vary_capital==1) then
        rebate_level=0.0D0
        lambdak=parameters(rebate_loop,1)
        rebate_slope=parameters(rebate_loop,2)
    elseif(vary_capital==2) then
        rebate_slope=0.0D0
        lambdak=parameters(rebate_loop,2)
        rebate_level=parameters(rebate_loop,1)
    else
        rebate_level=parameters(rebate_loop,1)
        rebate_slope=parameters(rebate_loop,2)
    end if
end if
lambda1=parameters(rebate_loop,3)+lambda1_benchmark

DO WHILE (diff_parameter>tol)
tax_diff=1.0D0
    DO WHILE (tax_diff >tol)
            CALL MPI_Bcast( prices, 11, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            tr=prices(1)
            lambda0=prices(2)
            w=prices(3)
            r=prices(4)
            tauss=prices(5)
            if(routine==5) then
                rebate_level=prices(6)
            end if
!            ss=prices(6)
            ss_cap=prices(7)
            lump_sum=prices(8)
            avg_earnings=prices(9)
            lambdak=prices(10)
            pe=prices(11)



!Solve Working
        do s=Jnr-1,2,-1
            CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            call MPI_RECV( end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            allocate(labor_alloc(typesN,shocksN,end_row-start_row+1))
            allocate(kprime_alloc(typesN,shocksN,end_row-start_row+1))
            allocate(v_alloc(typesN,shocksN,end_row-start_row+1))
            labor_alloc=0.0D0
            kprime_alloc=0
            v_alloc=0.0D0

            do kit=start_row,end_row
                if(kgrid(kit)<0.0D0) then
                    rate_use=r/survive(s-1)
                else
                    rate_use=r
                end if
                do e_it=1,shocksN
                    do i=1,typesN
                        c_upper=c_upper_fnc(epsilon1(s,1),kgrid(kit),kgrid,kgridN,tr,&
                            w*shocks(e_it)*types(i),rate_use)
                        temp=int(MINLOC(c_upper,1,mask=c_upper.gt.0.015D0+ebar))
                        if (c_upper(kgridN)>0.015D0+ebar) then
                            temp = kgridN
                        end if
                        if (c_upper(2)<=0.015D0+ebar) then
                            kprime_alloc(i,e_it,kit-start_row+1)=0
                            v_alloc(i,e_it,kit-start_row+1)=-1.0D0/0.0D0
                        else
                            CALL working_valuer(s,kgrid(kit),temp,tr,rate_use,&
                            survive(s),&
                            w*shocks(e_it)*types(i)*epsilon1(s,1),shock_pi(e_it,:,s,i),&
                            v_in(i,:,1:temp),&
                            v_alloc(i,e_it,kit-start_row+1),kprime_alloc(i,e_it,kit-start_row+1),&
                            labor_alloc(i,e_it,kit-start_row+1),&
                            partprod,hbar)


                        end if
                    end do
                end do
            end do
            call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( labor_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( kprime_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( v_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            deallocate(labor_alloc)
            deallocate(kprime_alloc)
            deallocate(v_alloc)
        end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Simulating
    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    call MPI_RECV(end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    allocate(sim_labor_alloc(J,end_row-start_row+1))
    allocate(sim_earnings_alloc(J,end_row-start_row+1))
    allocate(sim_capital_alloc(J+1,end_row-start_row+1))
    allocate(sim_consumption_alloc(J,end_row-start_row+1))
    allocate(sim_energy_alloc(J,end_row-start_row+1))
    allocate(sim_v_alloc(end_row-start_row+1))
    allocate(sim_total_income_alloc(J,end_row-start_row+1))
    allocate(sim_after_tax_income_alloc(J,end_row-start_row+1))
    !!Initialize

    sim_capital_alloc(1,:)=0.0D0
    sim_labor_alloc=0.0D0
    sim_earnings_alloc=0.0D0
    sim_total_income_alloc=0.0D0
    sim_after_tax_income_alloc=0.0D0




    do it=start_row,end_row
        !Working
        do s=1,Jnr-1
            if(sim_capital_alloc(s,it-start_row+1)<0.0D0) then
                if(s==1) then
                    rate_use=r/survive(1)
                else
                    rate_use=r/survive(s-1)
                end if
            else
                rate_use=r
            end if
            CALL working_valuer(s,sim_capital_alloc(s,it-start_row+1),kgridN,tr,rate_use,survive(s),&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1),&
                    shock_pi(sim_shocks(s,it),:,s,sim_types(s,it)),&
                    v(sim_types(s,it),:,:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),&
                    sim_labor_alloc(s,it-start_row+1),partprod,hbar)
            if (s==1) then
                sim_v_alloc(it-start_row+1)=dummy
            end if

            if (sim_labor_alloc(s,it-start_row+1)<hbar) then
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod)
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/(pe+taue)
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor_alloc(s,it-start_row+1)
            else
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/(pe+taue)
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor_alloc(s,it-start_row+1)
            end if
            sim_total_income_alloc(s,it-start_row+1)=sim_earnings_alloc(s,it-start_row+1)+&
            (sim_capital_alloc(s,it-start_row+1)+tr)*rate_use
            tax_dummy=all_taxer((sim_capital_alloc(s,it-start_row+1)+tr)*rate_use,sim_earnings_alloc(s,it-start_row+1),1)+taue*sim_energy_alloc(s,it-start_row+1)
            sim_after_tax_income_alloc(s,it-start_row+1)=sim_total_income_alloc(s,it-start_row+1)-tax_dummy
        end do



        do s=Jnr,J
            sim_labor_alloc(s,it-start_row+1)=0.0D0
            con=-1.0D0
            if(sim_shocks(Jnr,it) .le. shocksN/2) then
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
            else
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
            end if

            dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
            sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/(pe+taue)
            sim_total_income_alloc(s,it-start_row+1)=0.0D0+&
            (sim_capital_alloc(s,it-start_row+1)+tr)*r
            tax_dummy=all_taxer((sim_capital_alloc(s,it-start_row+1)+tr)*r,0.0D0,0)+taue*sim_energy_alloc(s,it-start_row+1)
            sim_after_tax_income_alloc(s,it-start_row+1)=sim_total_income_alloc(s,it-start_row+1)
        end do
    end do

    call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_earnings_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_alloc,(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_consumption_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_energy_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_v_alloc,(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_total_income_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_after_tax_income_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)


    deallocate(sim_labor_alloc)
    deallocate(sim_capital_alloc)
    deallocate(sim_energy_alloc)
    deallocate(sim_earnings_alloc)
    deallocate(sim_consumption_alloc)
    deallocate(sim_v_alloc)
    deallocate(sim_after_tax_income_alloc)
    deallocate(sim_total_income_alloc)

    CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end do
call MPI_Bcast( diff_parameter,1,MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
end do

end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Slaves that are not going to work for some
else if (my_id > (num_procs_run-1)) then
do rebate_loop=1,slopeN*levelN*progressivityN+(1+progressivityN)*extra
if(routine==5) then
    rebate_level=rebate_level_ig
else
    if(vary_capital==1) then
        rebate_level=0.0D0
        lambdak=parameters(rebate_loop,1)
        rebate_slope=parameters(rebate_loop,2)
    elseif(vary_capital==2) then
        rebate_slope=0.0D0
        lambdak=parameters(rebate_loop,2)
        rebate_level=parameters(rebate_loop,1)
    else
        rebate_level=parameters(rebate_loop,1)
        rebate_slope=parameters(rebate_loop,2)
    end if
end if
lambda1=parameters(rebate_loop,3)+lambda1_benchmark

!!!!!!!!!!!!!!!!!!!!!!!! Initialize loop stuff
out_loop=1
tax_loop = 1
tax_diff=1
non_clear=0
bad_combination_inner_stop=0
bad_combination_inner_nostop=0
bad_combination_outer=0
diff_parameter=1.0D0


DO WHILE (diff_parameter>tol)
tax_diff=1.0D0
    DO WHILE (tax_diff >tol)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Solving Value Function (These are not working)

            CALL MPI_Bcast( prices, 11, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            CALL MPI_Bcast( ss, typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            tr=prices(1)
            lambda0=prices(2)
            w=prices(3)
            r=prices(4)
            tauss=prices(5)
            if(routine==5) then
                rebate_level=prices(6)
            end if
!            ss=prices(6)
            ss_cap=prices(7)
            lump_sum=prices(8)
            avg_earnings=prices(9)
            lambdak=prices(10)
            pe=prices(11)


        do s=Jnr-1,2,-1
            CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
       end do




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Simulating
    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    call MPI_RECV(end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    allocate(sim_labor_alloc(J,end_row-start_row+1))
    allocate(sim_earnings_alloc(J,end_row-start_row+1))
    allocate(sim_capital_alloc(J+1,end_row-start_row+1))
    allocate(sim_consumption_alloc(J,end_row-start_row+1))
    allocate(sim_energy_alloc(J,end_row-start_row+1))
    allocate(sim_v_alloc(end_row-start_row+1))
    allocate(sim_total_income_alloc(J,end_row-start_row+1))
    allocate(sim_after_tax_income_alloc(J,end_row-start_row+1))
    !!Initialize

    sim_capital_alloc(1,:)=0.0D0
    sim_labor_alloc=0.0D0
    sim_earnings_alloc=0.0D0
    sim_total_income_alloc=0.0D0
    sim_after_tax_income_alloc=0.0D0




    do it=start_row,end_row
        !Working
        do s=1,Jnr-1
            if(sim_capital_alloc(s,it-start_row+1)<0.0D0) then
                if(s==1) then
                    rate_use=r/survive(1)
                else
                    rate_use=r/survive(s-1)
                end if
            else
                rate_use=r
            end if
            CALL working_valuer(s,sim_capital_alloc(s,it-start_row+1),kgridN,tr,rate_use,survive(s),&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1),&
                    shock_pi(sim_shocks(s,it),:,s,sim_types(s,it)),&
                    v(sim_types(s,it),:,:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),&
                    sim_labor_alloc(s,it-start_row+1),partprod,hbar)
            if (s==1) then
                sim_v_alloc(it-start_row+1)=dummy
            end if

            if (sim_labor_alloc(s,it-start_row+1)<hbar) then
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod)
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/(pe+taue)
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*partprod*sim_labor_alloc(s,it-start_row+1)
            else
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s,it-start_row+1),sim_capital_alloc(s+1,it-start_row+1),tr,rate_use,&
                    w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/(pe+taue)
                sim_earnings_alloc(s,it-start_row+1)=w*shocks(sim_shocks(s,it))*types(sim_types(s,it))*epsilon1(s,1)*sim_labor_alloc(s,it-start_row+1)
            end if
            sim_total_income_alloc(s,it-start_row+1)=sim_earnings_alloc(s,it-start_row+1)+&
            (sim_capital_alloc(s,it-start_row+1)+tr)*rate_use
            tax_dummy=all_taxer((sim_capital_alloc(s,it-start_row+1)+tr)*rate_use,sim_earnings_alloc(s,it-start_row+1),1)+taue*sim_energy_alloc(s,it-start_row+1)
            sim_after_tax_income_alloc(s,it-start_row+1)=sim_total_income_alloc(s,it-start_row+1)-tax_dummy
        end do



        do s=Jnr,J
            sim_labor_alloc(s,it-start_row+1)=0.0D0
            con=-1.0D0
            if(sim_shocks(Jnr,it) .le. shocksN/2) then
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)),avg_earnings)
            else
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1),&
                    kgrid(capital_positive:),&
                    tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),tr,r,survive(s),&
                    v_retired(sim_types(s,it),sim_shocks(s,it),:,s+1),dummy,sim_capital_alloc(s+1,it-start_row+1),avg_earnings,s)
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1),&
                    sim_capital_alloc(s+1,it-start_row+1),tr,r,ss(sim_types(s,it),sim_shocks(s,it)-shocksN/2),avg_earnings)
            end if

            dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1),con_dummy,s)
            sim_energy_alloc(s,it-start_row+1)=(con_dummy-sim_consumption_alloc(s,it-start_row+1))/(pe+taue)
            sim_total_income_alloc(s,it-start_row+1)=0.0D0+&
            (sim_capital_alloc(s,it-start_row+1)+tr)*r
            tax_dummy=all_taxer((sim_capital_alloc(s,it-start_row+1)+tr)*r,0.0D0,0)+taue*sim_energy_alloc(s,it-start_row+1)
            sim_after_tax_income_alloc(s,it-start_row+1)=sim_total_income_alloc(s,it-start_row+1)-tax_dummy
!            +&
!                rebate_calculator(0.0D0,(sim_capital_alloc(s,it-start_row+1)+tr)*r,0,0.0D0,tax_dummy)
        end do
    end do

    call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_earnings_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_alloc,(J+1)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_consumption_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_energy_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_v_alloc,(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_total_income_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_after_tax_income_alloc,(J)*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)


    deallocate(sim_labor_alloc)
    deallocate(sim_capital_alloc)
    deallocate(sim_energy_alloc)
    deallocate(sim_earnings_alloc)
    deallocate(sim_consumption_alloc)
    deallocate(sim_v_alloc)
    deallocate(sim_after_tax_income_alloc)
    deallocate(sim_total_income_alloc)

    CALL MPI_Bcast( tax_diff, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end do
call MPI_Bcast( diff_parameter,1,MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
end do
end do





end if
deallocate(population)








call MPI_FINALIZE(ierr)

end PROGRAM hetero_ss
